********************************************************************
****************** THE ILLUSION OF STABLE FERTILITY PREFERENCES ****
********************************************************************

clear all

global folder `c(pwd)'
global data "$folder\data"
global output "$folder\output_mainbody"
global logfiles "$folder\logfiles"
global fulldata "$data\ReplicationData_IofSFP_PS.dta"

* 1st time: ssc install labutil
set more off
 
******************************************************
******************* SETTINGS *************************
******************************************************
log using "$logfiles\logfile_mainbody", replace


*************************************************
*************** TABLE 1 *************************
** SURVEY TIMING & DATA AVAILABILITY ************
*************************************************
use "$fulldata", clear

**** ANALYSIS SAMPLE: ****
tab women_KLPS1 // number respondents Round 1
sum age_KLPS1 if women_KLPS1==1, detail // median age Round 1

tab women_KLPS12 // number respondents round 2
sum age_KLPS2 if women_KLPS12==1, detail // median age Round 2

tab women_KLPS13 // numbe respondents round 3 
sum age_KLPS3 if women_KLPS13==1, detail // median age round 3


**** EXTENDED SAMPLE ****
tab in_KLPS1 if female==1 // number respondents Round 1
sum age_KLPS1 if in_KLPS1==1 & female==1, detail // median age Round 1

tab in_KLPS2 if female==1 // number respondents Round 2
sum age_KLPS2 if in_KLPS2==1 & female==1, detail // median age Round 2

tab in_KLPS3 if female==1 // number respondents Round 3
sum age_KLPS3 if in_KLPS3==1 & female==1, detail // median age Round 3



*************************************************
*************** TABLE 2 *************************
** SUMMARY STATISTICS ***************************
*************************************************

**** ANALYSIS SAMPLE ****
use "$fulldata", clear
keep if women_KLPS1to3==1 // women in analysis sample (KLPS-1, 2, & 3)
rename weight_adj_KLPS1F weight_adj_KLPS1

forvalues i=1(1)3 {
preserve
foreach var in age intchi living_kids ind_parent married { 
	rename `var'_KLPS`i' `var'_summ
	}
eststo women_KLPS`i'F: estpost sum age_summ intchi_summ living_kids_summ ind_parent_summ married_summ [aw=weight_adj_KLPS`i']
restore 
}

/* for KLPS-1: for each variable, there is at least one respondent missing data --> N=238. But N=239 participated in all 3 rounds, and for 
N=239 there is information for at least one of the variables */ 


**** EXTENDED SAMPLE ****
use "$fulldata", clear
gen weight_adj_KLPS1=weight_adj_KLPS1H
forvalues i=1(1)3 {
preserve
foreach var in age intchi living_kids ind_parent married { 
	rename `var'_KLPS`i' `var'_summ
	}
eststo women_KLPS`i': estpost summarize age_summ intchi_summ living_kids_summ ind_parent_summ married_summ [aw=weight_adj_KLPS`i'] if female==1
restore
}

rename age_KLPS1 age_summ
rename intchi_KLPS1 intchi_summ 
rename living_kids_KLPS1 living_kids_summ
rename ind_parent_KLPS1 ind_parent_summ
rename married_KLPS1 married_summ

esttab women_KLPS1F women_KLPS2F women_KLPS3F women_KLPS1 women_KLPS2 women_KLPS3 ///
using "$output/table2.tex", replace f nomtitle label collabels(none) cells("mean(pattern(1 1 1) fmt(2))")


**************************************************************************************
*************** TABLE 3 **************************************************************
** REGRESSIONS OF ACTUAL FERTILITY ON REPRODUCTIVE DESIRES ***************************
**************************************************************************************
use "$fulldata", clear
keep if women_KLPS1to3==1 // women in analysis sample (KLPS-1, 2, & 3)

*** ALL WOMEN ***
forvalues i=2(1)3 {
sum live_births_K1K`i' [aw=weight_adj_KLPS`i'] if addintchi_KLPS1!=. // Weighted mean # of additional children between rounds (+ std. dev.)
reg live_births_K1K`i' addintchi_KLPS1 [pw=weight_adj_KLPS`i'], vce(cluster base_schid) // use intensive tracking adjusted weights and cluster at baseline school level
eststo reg`i'
}

*** Pregnancies>0 (Rd 1) ***
forvalues i=2(1)3 {
sum live_births_K1K`i' [aw=weight_adj_KLPS`i'] if ind_pregnancy_KLPS1==1 & addintchi_KLPS1!=. // Weighted mean # of additional children between rounds
reg live_births_K1K`i' addintchi_KLPS1 [pw=weight_adj_KLPS`i'] if ind_pregnancy_KLPS1==1, vce(cluster base_schid)
eststo regadd`i'
}

*** Never Pregnant (Rd 1) ***
forvalues i=2(1)3 {
sum live_births_KLPS`i' [aw=weight_adj_KLPS`i'] if ind_pregnancy_KLPS1==0 & addintchi_KLPS1!=. // Weighted mean # of additional children between rounds
reg live_births_KLPS`i' addintchi_KLPS1 [pw=weight_adj_KLPS`i'] if ind_pregnancy_KLPS1==0, vce(cluster base_schid)
eststo reg`i'_nonpregnant
}


*** Gather coefficients, N, & R-squared to embed in table
esttab reg2 reg3 regadd2 regadd3 reg2_nonpregnant reg3_nonpregnant ///
	using "$output\table3.tex", ///
	replace f collabels(none) b(3) se(3) star(* 0.10 ** 0.05 *** 0.01) label ///
	nonumbers nomtitle ///
	title("Dep. Var.: Additional Children between Round 1 and ...") keep(addintchi_KLPS1) ///
	se stats(N r2, fmt(%9.0g %9.3f) labels(N R-squared)) 	


************************************************************************************************
*************** FIGURE 1 ***********************************************************************
** DISTRIBUTION OF CHANGES IN DESIRED CHILDREN BETWEEN SURVEY ROUNDS ***************************
************************************************************************************************
use "$fulldata", clear
keep if women_KLPS1to3==1 // women in analysis sample (KLPS-1, 2, & 3)
foreach var in D_pref_KLPS2_1 D_pref_KLPS3_2 D_pref_KLPS3_1 {
keep if `var'>=-4 & `var'<=4
}

* Histogram of changes between survey rounds 1 and 2
preserve
gen weight_rounded=round(weight_adj_KLPS2) // use intensive survey tracking adjusted weights of later survey round
expand weight_rounded
graph set window fontface "Times New Roman"
sum D_pref_KLPS2_1
local avg1 : di %03.2f `r(mean)'
histogram D_pref_KLPS2_1, discrete width(1) ///
	xlabel (-4 "-4" -2 "-2" 0 "0" 2 "2" 4 "4") 	xtitle("Change in desired number of children", size(small)) xsc(titlegap(1.5)) ///
	ylabel (0(0.1)0.5, format(%03.1f)) ytitle("Density", size(small)) ///
	title("Round 1 to 2" "(ca age 19 to 23)", color(black)) name(D_pref_KLPS2_1_w, replace) ///
	fcolor(dknavy) lcolor(white) graphregion(color(white)) xline(`r(mean)') ///
	text(0.45 2.4 "Mean: `avg1'", color(red)) 
restore

* Histogram of changes between survey rounds 2 and 3
preserve
gen weight_rounded=round(weight_adj_KLPS3) // use intensive survey tracking adjusted weights of later survey round
expand weight_rounded

sum D_pref_KLPS3_2
local avg1 : di %03.2f `r(mean)'
histogram D_pref_KLPS3_2, discrete width(1) ///
	xlabel(-4 "-4" -2 "-2" 0 "0" 2 "2" 4 "4") ///
	xtitle("Change in desired number of children", size(small)) xsc(titlegap(1.5)) ///
	ylabel (0(0.1)0.5, format(%03.1f)) ytitle("Density", size(small)) ///
	title("Round 2 to 3" "(ca age 23 to 28)", color(black)) name(D_pref_KLPS3_2_w, replace) ///
	fcolor(dknavy) lcolor(white) graphregion(color(white)) xline(`r(mean)') ///
	text(0.45 2.4 "Mean: `avg1'", color(red))

* Histogram of changes between survey rounds 1 and 3
sum D_pref_KLPS3_1
local avg1 : di %03.2f `r(mean)'
histogram D_pref_KLPS3_1, discrete width(1) ///
	xlabel(-4 "-4" -2 "-2" 0 "0" 2 "2" 4 "4") ///
	xtitle("Change in desired number of children", size(small)) xsc(titlegap(1.5)) ///
	ylabel (0(0.1)0.5, format(%03.1f)) ytitle("Density", size(small)) ///
	title("Round 1 to 3" "(ca age 19 to 28)", color(black)) name(D_pref_KLPS3_1_w, replace) ///
	fcolor(dknavy) lcolor(white) graphregion(color(white)) xline(`r(mean)') ///
	text(0.45 2.4 "Mean: `avg1'", color(red))

* Combine histograms to final graph
graph combine D_pref_KLPS2_1_w D_pref_KLPS3_2_w D_pref_KLPS3_1_w, rows(1) ycommon xcommon ///
title("Change in the desired number of children", color(black)) graphregion(color(white))

graph export "$output\figure1.pdf", as(pdf) replace
restore


**********************************************************************************************
*************** FIGURE 2 *********************************************************************
** DISTRIBUTION OF DESIRED NUMBER OF CHILDREN ACROSS SURVEY ROUNDS ***************************
**********************************************************************************************
graph set window fontface "Times New Roman"

use "$fulldata", clear
keep if women_KLPS1to3==1  // women in analysis sample (KLPS-1, 2, & 3)
keep if intchi_KLPS1>=0 & intchi_KLPS1<=8 & intchi_KLPS2<=8 & intchi_KLPS2>=0

gen on_diagonal_D_KLPS2_1 = 0 
replace on_diagonal_D_KLPS2_1 = 1 if intchi_KLPS2==intchi_KLPS1
gen above_diagonal_D_KLPS2_1 = 0 
replace above_diagonal_D_KLPS2_1 = 1 if intchi_KLPS2>intchi_KLPS1
gen below_diagonal_D_KLPS2_1 = 0 
replace below_diagonal_D_KLPS2_1 = 1 if intchi_KLPS2<intchi_KLPS1

sum on_diagonal_D_KLPS2_1 // unweighted (Number of observations)
local N = `r(N)'

gen weight_rounded=round(weight_adj_KLPS2)
expand weight_rounded

sum on_diagonal_D_KLPS2_1 // weighted
local ondiagonal : di %3.1f `r(mean)'*100
local offdiagonal = 100-`ondiagonal'
sum above_diagonal_D_KLPS2_1
local abovediagonal : di %3.1f `r(mean)'*100
sum below_diagonal_D_KLPS2_1 
local belowdiagonal : di %3.1f `r(mean)'*100

egen obs_changes_KLPS12 = count(pupid), by(intchi_KLPS1 intchi_KLPS2)

twoway (scatter intchi_KLPS2 intchi_KLPS1 [w=obs_changes_KLPS12], msymbol(circle_hollow) ytitle("Desired number of children at Round 2") ///
xtitle("Desired number of children at Round 1") xlabel(0(1)8, valuelabels) ylabel(0(1)8, valuelabels angle(0)) graphregion(color(white))) ///
(scatteri 0 0 6.5 6.5, recast(line) legend(off)) (function y=7.7, range(6.5 9.7) color(black)) ///
, name(KLPS12, replace) ysc(range(0 8)) xsc(range(0 8)) ///
text(8 7.2 "Sample size", just(left)  margin(l+2 t+1 b+1))  ///
text(7.3 7.3 "On-diagonal" "Off-diagonal", just(left)  margin(l+2 t+1 b+1))  ///
text(6.6 7.6 "{it:Below diagonal}" "{it:Above diagonal}",  just(left) margin(l+2 t+1 b+1)) ///
text(8 9.2 "`N'", just(right)  margin(l+2 t+1 b+1))  ///
text(7.3 9.2 "`ondiagonal'%" "`offdiagonal'%", just(right) margin(l+2 t+1 b+1)) ///
text(6.6 9.2 "{it:`belowdiagonal'%}" "{it:`abovediagonal'%}", just(right) margin(l+2 t+1 b+1))

graph export "$output\figure2a.pdf", as(pdf) replace

* CHANGES KLPS-2 to KLPS-3
use "$fulldata", clear
keep if women_KLPS1to3==1 & intchi_KLPS2>=0 & intchi_KLPS2<=8 & intchi_KLPS3<=8 & intchi_KLPS3>=0

gen on_diagonal_D_KLPS3_2 = 0 
replace on_diagonal_D_KLPS3_2 = 1 if intchi_KLPS3==intchi_KLPS2
gen above_diagonal_D_KLPS3_2 = 0 
replace above_diagonal_D_KLPS3_2 = 1 if intchi_KLPS3>intchi_KLPS2
gen below_diagonal_D_KLPS3_2 = 0 
replace below_diagonal_D_KLPS3_2 = 1 if intchi_KLPS3<intchi_KLPS2

sum on_diagonal_D_KLPS3_2 // unweighted (Number of observations)
local N = `r(N)'

gen weight_rounded=round(weight_adj_KLPS3)
expand weight_rounded

sum on_diagonal_D_KLPS3_2 // weighted
local ondiagonal : di %3.1f `r(mean)'*100
local offdiagonal = 100-`ondiagonal'
sum above_diagonal_D_KLPS3_2 
local abovediagonal : di %3.1f `r(mean)'*100
sum below_diagonal_D_KLPS3_2
local belowdiagonal : di %3.1f `r(mean)'*100

egen obs_changes_KLPS23 = count(pupid), by(intchi_KLPS2 intchi_KLPS3)

twoway (scatter intchi_KLPS3 intchi_KLPS2 [w=obs_changes_KLPS23], msymbol(circle_hollow) ytitle("Desired number of children at Round 3") ///
xtitle("Desired number of children at Round 2") xlabel(0(1)8, valuelabels) ylabel(0(1)8, valuelabels angle(0)) graphregion(color(white))) ///
(scatteri 0 0 6.5 6.5, recast(line) legend(off)) (function y=7.7, range(6.5 9.7) color(black)) ///
, name(KLPS12, replace) ysc(range(0 8)) xsc(range(0 8)) ///
text(8 7.2 "Sample size", just(left)  margin(l+2 t+1 b+1))  ///
text(7.3 7.3 "On-diagonal" "Off-diagonal", just(left)  margin(l+2 t+1 b+1))  ///
text(6.6 7.6 "{it:Below diagonal}" "{it:Above diagonal}",  just(left) margin(l+2 t+1 b+1)) ///
text(8 9.2 "`N'", just(right)  margin(l+2 t+1 b+1))  ///
text(7.3 9.2 "`ondiagonal'%" "`offdiagonal'%", just(right) margin(l+2 t+1 b+1)) ///
text(6.6 9.2 "{it:`belowdiagonal'%}" "{it:`abovediagonal'%}", just(right) margin(l+2 t+1 b+1))

graph export "$output\figure2b.pdf", as(pdf) replace

* CHANGES KLPS-1 to KLPS-3
use "$fulldata", clear
keep if women_KLPS1to3==1 & intchi_KLPS1>=0 & intchi_KLPS1<=8 & intchi_KLPS3<=8 & intchi_KLPS3>=0

gen on_diagonal_D_KLPS3_1 = 0 
replace on_diagonal_D_KLPS3_1 = 1 if intchi_KLPS3==intchi_KLPS1
gen above_diagonal_D_KLPS3_1 = 0 
replace above_diagonal_D_KLPS3_1 = 1 if intchi_KLPS3>intchi_KLPS1
gen below_diagonal_D_KLPS3_1 = 0 
replace below_diagonal_D_KLPS3_1 = 1 if intchi_KLPS3<intchi_KLPS1

sum on_diagonal_D_KLPS3_1 // unweighted (Number of observations)
local N = `r(N)'

gen weight_rounded=round(weight_adj_KLPS3)
expand weight_rounded

sum on_diagonal_D_KLPS3_1 // weighted
local ondiagonal : di %3.1f `r(mean)'*100
local offdiagonal = 100-`ondiagonal'
sum above_diagonal_D_KLPS3_1 
local abovediagonal : di %3.1f `r(mean)'*100
sum below_diagonal_D_KLPS3_1
local belowdiagonal : di %3.1f `r(mean)'*100

egen obs_changes_KLPS13 = count(pupid), by(intchi_KLPS1 intchi_KLPS3)

twoway (scatter intchi_KLPS3 intchi_KLPS1 [w=obs_changes_KLPS13], msymbol(circle_hollow) ytitle("Desired number of children at Round 3") ///
xtitle("Desired number of children at Round 1") xlabel(0(1)8, valuelabels) ylabel(0(1)8, valuelabels angle(0)) graphregion(color(white))) ///
(scatteri 0 0 6.5 6.5, recast(line) legend(off)) (function y=7.7, range(6.5 9.7) color(black)) ///
, name(KLPS12, replace) ysc(range(0 8)) xsc(range(0 8)) ///
text(8 7.2 "Sample size", just(left)  margin(l+2 t+1 b+1))  ///
text(7.3 7.3 "On-diagonal" "Off-diagonal", just(left)  margin(l+2 t+1 b+1))  ///
text(6.6 7.6 "{it:Below diagonal}" "{it:Above diagonal}",  just(left) margin(l+2 t+1 b+1)) ///
text(8 9.2 "`N'", just(right)  margin(l+2 t+1 b+1))  ///
text(7.3 9.2 "`ondiagonal'%" "`offdiagonal'%", just(right) margin(l+2 t+1 b+1)) ///
text(6.6 9.2 "{it:`belowdiagonal'%}" "{it:`abovediagonal'%}", just(right) margin(l+2 t+1 b+1))

graph export "$output\figure2c.pdf", as(pdf) replace


*******************************************************************
*************** FIGURE 3 ******************************************
** EXPECTATIONS FOR DIFFERENT SCENARIOS ***************************
*******************************************************************
use "$fulldata", clear

keep if women_KLPS1==1 // Women in analysis-sample (KLPS-1)
gen weight_rounded=round(weight_adj_KLPS1F)
expand weight_rounded
#delimit;
local scenarios "k1f_s5_14b_marrysoon k1f_s5_14c_nohusband k1f_s5_14d_allfemale 
k1f_s5_14e_lessincome k1f_s5_14f_death k1f_s5_14g_givenanotherchild k1f_s5_14h_daughterdelivers 
k1f_s5_14i_difficultpreg k1f_s5_14j_allmale k1f_s5_14k_disagreement k1f_s5_14l_increasedincome 
k1f_s5_14m_cowifemanychildren k1f_s5_14n_3morekids k1f_s5_14o_daughterinlawdelivers 
k1f_s5_14p_manneedsmore k1f_s5_14q_cowifeleaves k1f_s5_14r_husbandmarries k1f_s5_14s_1kidgrowsoutsidehome";
#delimit cr // left out: k1f_s5_14a_otherwives
local i=1
tab k1f_s5_14a_otherwives, matcell(x`i')
local name: variable label k1f_s5_14a_otherwives
matrix list x`i'
clear 
set obs 1
gen question = "`name'" in 1
gen code = 1 in 1
gen answer_same = x`i'[1,1]
gen answer_more =x`i'[2,1]
gen answer_fewer = x`i'[3,1]
gen answer_NA = x`i'[4,1]
gen answer_DK = x`i'[5,1]
save "$output\scenarios.dta", replace

foreach var in `scenarios' {
use "$fulldata", clear
keep if women_KLPS1==1
gen weight_rounded=round(weight_adj_KLPS1F)
expand weight_rounded
local name : variable label `var' 
local i=`i'+1
tab `var', matcell(x`i')
use "$output\scenarios.dta", clear
set obs `i'
replace question = "`name'" in `i'
replace code = `i' in `i'
replace answer_same = x`i'[1,1] in `i'
replace answer_more =x`i'[2,1] in `i'
replace answer_fewer = x`i'[3,1] in `i'
replace answer_NA = x`i'[4,1] in `i'
replace answer_DK = x`i'[5,1] in `i'
save "$output\scenarios.dta", replace
sleep 200
}

use "$output\scenarios.dta", clear
egen total_answers = rowtotal(answer_same answer_more answer_fewer)
foreach var in answer_same answer_more answer_fewer {
replace `var' = `var'/total_answers
}
sort answer_same

* ORDER: FINANCES & HEALTH
gen label_order = 1 if code == 12 // finances improve
replace label_order = 2 if code == 5 // finances worsen
replace label_order = 3 if code == 9 // pregnancies are difficult

* ORDER: MARRIAGES & HOUSEHOLD
replace label_order = 4 if code == 16 // husband wants more children
replace label_order = 5 if code == 17 // left alone with husband (co-wife leaves)
replace label_order = 6 if code == 2 // marry soon
replace label_order = 7 if code == 18 // husband takes another wife
replace label_order = 8 if code == 13 // co-wife has many children
replace label_order = 9 if code == 1 // become a junior co-wife
replace label_order = 10 if code == 11 // no longer get along with spouse
replace label_order = 11 if code == 3 // unable to find husband

* ORDER: CHILDREN
replace label_order = 12 if code == 19 // child fostered away
replace label_order = 13 if code == 4 // all desired children female
replace label_order = 14 if code == 10 // all desired children male
replace label_order = 15 if code == 6 // a child dies in infancy
replace label_order = 16 if code == 7 // receive a teen foster child
replace label_order = 17 if code == 14 // receive 3 young foster children
replace label_order = 18 if code == 15 // daughter-in-law gives birth
replace label_order = 19 if code == 8 // daughter gives birth

labmask label_order, values(question)

graph hbar (mean) answer_more answer_same answer_fewer, over(label_order) ///
stack bar(1, color(gs7) lcolor(gs7)) bar(2, color(gs15) lcolor(gs15)) ///
bar(3, color(gs1) lcolor(gs1)) graphregion(color(white)) ///
legend(order(1 "More" 2 "Same" 3 "Less") rows(1)) 

graph export "$output\Figure3.pdf", as(pdf) replace


**************************************************
*************** FIGURE 4 *************************
***** RECALL PATTERNS ****************************
**************************************************
use "$fulldata", clear
keep if women_KLPS12==1 // women in analysis sample (KLPS-1 & KLPS-2)

* generate indicator of respondents recalling having lowered, increased, or kept stable desires
gen recallslowered = 0 if in_KLPS2==1 & intchi_KLPS2!=. & k2_s17_2_15oldideal!=.
replace recallslowered = 1 if k2_s17_2_15oldideal>intchi_KLPS2 & intchi_KLPS2!=. & k2_s17_2_15oldideal!=.
gen recallsincreased = 0 if in_KLPS2==1 & intchi_KLPS2!=. & k2_s17_2_15oldideal!=.
replace recallsincreased = 1 if k2_s17_2_15oldideal<intchi_KLPS2 & intchi_KLPS2!=. & k2_s17_2_15oldideal!=.
gen recallsstable = 0 if in_KLPS2==1 & intchi_KLPS2!=. & k2_s17_2_15oldideal!=.
replace recallsstable = 1 if k2_s17_2_15oldideal==intchi_KLPS2 & intchi_KLPS2!=. & k2_s17_2_15oldideal!=.

* generate variable indicating the direction of recalled change in desires
gen recalleddirection = -1 if recallslowered==1
replace recalleddirection = 0 if recallsstable==1
replace recalleddirection = 1 if recallsincreased==1

* generate variable indicating the actual direction of change in stated desires
gen direction_change = . if women_KLPS12==1
replace direction_change = 0 if D_pref_KLPS2_1<0 
replace direction_change = 1 if D_pref_KLPS2_1==0
replace direction_change = 2 if D_pref_KLPS2_1>0 & D_pref_KLPS2_1!=.
label define change2 0 "Lowered desires" 1 "Stable desires" 2 "Increased desires" 
label values direction_change change2

tab direction_change // unweighted distribution 

gen weight_rounded=round(weight_adj_KLPS2) // use intensive-tracking adjusted weights from later survey round (round 2)
expand weight_rounded

tab direction_change // weighted distribution

* Encode the share of respondents having lowered, increased, or kept stable desires
gen sharelowered = 0 if direction_change!=.
replace sharelowered = 1 if direction_change==0
sum sharelowered 
local sharelowered = `r(mean)'

gen shareincreased = 0 if direction_change!=.
replace shareincreased = 1 if direction_change==2
sum shareincreased
local shareincreased = `r(mean)'

gen sharestable = 0 if direction_change!=.
replace sharestable = 1 if direction_change==1
sum sharestable
local sharestable = `r(mean)'

* Encode the share of respondents recalling having lowered desires conditional on actual change in desires
sum recallslowered if direction_change==0
local shareRlowered_L = `r(mean)'

sum recallslowered if direction_change==1
local shareRlowered_S = `r(mean)'

sum recallslowered if direction_change==2
local shareRlowered_I = `r(mean)'

* Encode the share of respondents recalling having increased desires conditional on actual change in desires
sum recallsincreased if direction_change==0
local shareRincreased_L = `r(mean)'

sum recallsincreased if direction_change==1
local shareRincreased_S = `r(mean)'

sum recallsincreased if direction_change==2
local shareRincreased_I = `r(mean)'

* Encode the share of respondents recalling having kept stable desires conditional on actual change in desires
sum recallsstable if direction_change==0
local shareRstable_L = `r(mean)'

sum recallsstable if direction_change==1
local shareRstable_S = `r(mean)'

sum recallsstable if direction_change==2
local shareRstable_I = `r(mean)'

* generate variable with overall shares (in whole sample) of those recalling having lowered desires for each actual direction of change
gen absolute_shareRlower = .
replace absolute_shareRlower = `sharelowered'*`shareRlowered_L' if direction_change==0
replace absolute_shareRlower = `sharestable'*`shareRlowered_S' if direction_change==1
replace absolute_shareRlower = `shareincreased'*`shareRlowered_I' if direction_change==2

* generate variable with overall shares (in whole sample) of those recalling having increased desires for each actual direction of change
gen absolute_shareRhigher = .
replace absolute_shareRhigher = `sharelowered'*`shareRincreased_L' if direction_change==0
replace absolute_shareRhigher = `sharestable'*`shareRincreased_S' if direction_change==1
replace absolute_shareRhigher = `shareincreased'*`shareRincreased_I' if direction_change==2

* generate variable with overall shares (in whole sample) of those recalling having kept stable desires for each actual direction of change
gen absolute_shareRstable = .
replace absolute_shareRstable = `sharelowered'*`shareRstable_L' if direction_change==0
replace absolute_shareRstable = `sharestable'*`shareRstable_S' if direction_change==1
replace absolute_shareRstable = `shareincreased'*`shareRstable_I' if direction_change==2

foreach t in L S I {
local shareRlowered2_`t' = round(100*`shareRlowered_`t'',1.0)
local shareRstable2_`t' = round(100*`shareRstable_`t'',1.0)
local shareRincreased2_`t' = round(100*`shareRincreased_`t'',1.0)
}

* Produce graph indicating share of individuals who lowered, increased and kept stable desires as well as memory conditional on the actual direction of change
graph hbar (mean) absolute_shareRlower absolute_shareRstable absolute_shareRhigher, over(direction_change) ///
stack bar(1, color(gs12%66) lcolor(gs10)) bar(2, color(gs6%66) lcolor(gs4)) ///
bar(3, color(gs0%66) lcolor(gs0)) graphregion(color(white)) ///
legend(order(1 "Recalls having lowered desires" 2 "Recalls stable desires" ///
	3 "Recalls having increased desires" 4 "{bf:Bold print marks the correct recall direction}") position(12) rows(4)) ///
text(0.02 101 "{bf:`shareRlowered2_L'%}", color(black) size(medlarge)) text(0.18 101 "`shareRstable2_L'%", color(black)) ///
text(0.31 101 "`shareRincreased2_L'%", color(black)) text(0.02 65 "`shareRlowered2_S'%", color(black)) ///
text(0.18 65 "{bf:`shareRstable2_S'%}", color(black) size(medlarge)) text(0.33 65 "`shareRincreased2_S'%", color(black)) ///
text(0.012 29 "`shareRlowered2_I'%", color(black)) text(0.14 29 "`shareRstable2_I'%", color(black)) ///
text(0.265 29 "{bf:`shareRincreased2_I'%}", color(black) size(medlarge)) ///
ytitle("Proportions (analysis sample)", placement(3)) ysc(titlegap(2)) text(-0.07 110 "Change in desires") ///
ylabel(0(0.1)0.4, format(%03.1f))

graph export "$output\Figure4.pdf", as(pdf) replace 

log close 
